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We consider a generalization of the standard Metropolis algorithm acceptance/rejection decision 
rule and numerically explore its properties using auxiliary field quantum Monte Carlo. The gen- 
eralization involves a free parameter which, given a criterion for proposing attempted moves, can 
be used to tune the average acceptance rate in a particular way. Such tuning can also potentially 
change Monte Carlo autocorrelation times, and the combination of the changing acceptance rate 
and autocorrelation times raises the possibility of more efficient simulations. We explore these issues 
using primarily massively parallel quantum Monte Carlo runs of the "test case" two-dimensional 
Hubbard model, and discuss results and applications. 
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I. INTRODUCTION 

Since its introduction H, the "Metropolis algorithm" 
for Monte Carlo has become a widely used and power- 
ful numerical tool, both for classical and quantum sys- 
tems. The algorithm describes a Markov chain through 
states of the system with rules for first proposing a new 
state and then deciding whether or not to accept a move 
to the proposed state. In statistical mechanics simula- 
tions, the Metropolis rules by construction sample states 
according to their Boltzmann weights, allowing compu- 
tation of expectation values. However, this Boltzmann 
weight sampling can be satisfied by other ways of mak- 
ing the acceptance/rejection decision besides the conven- 
tional one. We consider here a class of such alternative 
acceptance/rejection decisions and explore whether they 
might be used to increase simulation efficiency. We focus 
in this paper on auxiliary-field fermion quantum Monte 
Carlo simulations |PHl0| though our approach and results 
may have more general application, particularly in other 
electron quantum Monte Carlo methods [p~T] [Tqj] . We use 
primarily massively parallel machines in the quantum 
Monte Carlo simulations, which has allowed us to more 
easily gather data than would otherwise be possible. 

One way to increase general Monte Carlo efficiency is 
to minimize the autocorrelation times r e . r e is qualita- 
tively defined for an observable as the number of Monte 
Carlo steps over which measurements of that observable 
remain correlated once the system has equilibrated (i.e., 
once the system has become independent of the initial 
configuration). N Monte Carlo steps then correspond to 
approximately iV; n d = N/r e independent measurements. 
A particular statistical error will require a given number 
iV; n d of uncorrelated data points. Hence, N = N^Te 
steps will be required for a given desired statistical er- 
ror, with reducing r e commensurately reducing Monte 
Carlo run time. A related consideration involves the 
"warm-up" or equilibration time t w , a measure of the 



average number of steps required to go from typically 
low-weight initial configurations to the high-weight re- 
gion of the phase space where measurements can then be 
taken. 

In auxiliary field and other determinantal quantum 
Monte Carlo methods, however, another factor is intro- 
duced with regards to efficiency: an accepted move is 
often so much more computationally intensive than a re- 
jected move that the time required for rejected moves can 
be neglected compared to the computationally dominat- 
ing acceptance calculations [p|-p|,|i7[ . Then, the dominant 
computational time required for a simulation of N steps 
is proportional rather to NA e = N- m ^A e T e , where A e is 
the average equilibrated move acceptance rate, so that 
one now wishes to minimize the quantity A e T e . Anal- 
ogously, the dominant required computational time for 
equilibration will be A w t w , where A w is an acceptance 
rate averaged during the warm-up process. Exploring the 
above efficiency issues is a major focus of the paper. One 
aspect of this involves the calculation of equilibration and 
autocorrelation times for the "test case" two-dimensional 
Hubbard model, which results, to the best of our knowl- 
edge, have not previously appeared. 

After the introduction, we discuss various Monte Carlo 
algorithmic considerations, including parallelization is- 
sues, a generalization of the conventional Metropolis ac- 
ceptance/rejection rule, and relevance to auxiliary field 
and other determinantal quantum Monte Carlo. We then 
explore properties of the new Metropolis generalization 
using primarily massively parallel auxiliary field quantum 
Monte Carlo simulations of the two-dimensional Hubbard 
model. Lastly, we discuss the numerical results and sum- 
marize. 



II. MONTE CARLO ALGORITHM 
CONSIDERATIONS 
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A. Parallelization Issues 

Currently, the fastest computers available are those 
with a massively parallel architecture. These computers 
consist of hundreds to thousands of individual processors 
linked so that they can communicate with each other. 
Efficient use of these machines is governed by keeping 
communication time sufficiently low that computational 
speed scales roughly linearly with the number of proces- 
sors. 

One approach for using such computers in Monte Carlo 
simulations is to distribute a single Monte Carlo "walk" 
over all the processors. This allows the computational 
work of each step in the walk to be spread out over all 
the processors, but it can require a high degree of com- 
munication. This approach can be feasible for classical 
simulations of large systems with short-range forces, with 
each processor "owning" part of the system. However, 
it has not been shown to be efficient for determinantal 
quantum Monte Carlo simulations of the type we will 
discuss, and in particular might be expected to be less 
efficient for two dimensional systems. 

A second approach is to let every processor perform 
its own independent Monte Carlo "walk". Each proces- 
sor then starts with a different initial configuration and 
can proceed independently of the other processors, and 
no communication is required until the final accumula- 
tion of results 

A potential drawback to this second approach, how- 
ever, is that each processor must independently move 
from its typically low-weight initial configuration to 
the high- weight region of the phase space ("equilibra- 
tion" ) before measurements can start to be taken |l9|,^0| . 
Hence, while the measurement process is itself parallel, 
with measurements from the different processors com- 
bining to reduce statistical error, equilibration is serial. 
More specifically, let Ni a d denote the number of inde- 
pendent measurements required for a desired statistical 
error, let t w denote the number of Monte Carlo steps re- 
quired for equilibration ( "warm-up" ) , let r e denote auto- 
correlation time in Monte Carlo steps, and let Np denote 
the number of processors in a parallel machine. Then, for 
a serial or vector machine, the required number of Monte 
Carlo steps N is given by 



N = T W + N ind T e . 



(1) 



Assuming a fluctuation-dissipation-type scenario, so that 
t w and r e have similar values, equilibration then plays a 
small role in the required computational time [ pl| . How- 
ever, for a parallel machine, the number of steps required 
is 



N 



N ind T e /N F 



(2) 



Hence, for a massively parallel machine with thousands 
of processors, equilibration provides a potential bottle- 
neck, which could be improved by reducing t w . It seems 



that such a bottleneck may actually be alleviated some- 
what by the "sign problem" , which can necessitate a very 
large value of N ind for reasonable statistical error. 



B. Metropolis, Symmetric, and "Generalized" 
Decision Rules 

In Monte Carlo simulations, new moves are proposed 
according to some procedure and a decision is then made 
whether to accept the proposed move or to reject the 
move and remain in the current state. The proposal and 
decision together usually satisfy "microscopic reversibil- 
ity", or "detailed balance" p2fl . The most commonly 
used proposal procedure stipulates that the probability 
of proposing a move to state j given that one is currently 
in state i is identical to the probability of proposing a 
move to state i given that one is in state j, though other 
procedures have also been used |^3|j24| . Although our 
analysis could be generalized, we will assume the above 
"symmetric" move proposal procedure throughout this 
paper in the case of detailed balance. 

Probably the most common rule for deciding whether 
or not to accept a proposed move is the "Metropolis de- 
cision" . Let the Monte Carlo sampling be over the Boltz- 
mann distribution, let E denote the energy of the current 
state, let E' denote the energy of the proposed state, let 
AE = E' - E, and let /3 = l/(k B T), where T is the 
temperature and ks is Boltzmann's constant. Then, as 
in the original Metropolis paper [Q , the probability P of 
accepting the proposed state is given by 



P 



M 



,-0AE 
1 



, AE > 
AE < 0. 



(3) 



Another decision rule which has also been used is the 
so-called "symmetric rule" [p5| , p6| , 

e -/3AB 

P S = 0AE ■ ( 4 ) 



1 



e 



We note that, since Pm > Ps for all (3AE, Pm will give 
a higher average acceptance rate than P$- 

It was shown by Peskun under quite general assump- 
tions that the Metropolis decision would lead to smaller 
statistical errors in the limit of very long simulations 
than would the symmetric decision if detailed balance 
were satisfied p^ , p8| ], this result being most relevant to 
T e . The result correlates with the higher Metropolis de- 
cision acceptance rate. When looking at the convergence 
of state distributions "operated on" by Monte Carlo de- 
cision rules, however, it was found that the symmetric 
decision rule was superior in certain cases, particularly 
those where a small number of states was accessible at 
each Monte Carlo step and where the quantity j3AE was 
typically around mag nitude 1.0 or less ||S This lat " 
ter result is more directly relevant to t w . Further, the 
Peskun analysis does not apply to certain cases of inter- 
est which can lead to the correct limiting (e.g., Boltz- 
mann) distribution but which do not necessarily satisfy 
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detailed balance, such as systematically moving through 
the lattice of a discrete system when proposing moves as 
opposed to randomly selecting lattice sites Q . 

We explore here a gen eralization of the Metropolis and 
symmetric decisions |pl| , 



P = 



e -0AE 

r+ e -' 3AE 



AE > 
AE < 



(5) 



where T is a tunable parameter. Like the standard 
Metropolis and symmetric decision rules, the P of Eq. |^ 
satisfies the condition 

e -PE p( E ^ £/) = e -0E' p ( E > ^ E j (g) 

When r = we recover the Metropolis decision rule and 
when r = 1 we recover the symmetric rule. For < 
r < 1, P smoothly interpolates between the Metropolis 
and symmetric limits. However, P is also defined for any 
r > 1. 
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FIG. 1. Effect of F on acceptance probability P of Eq. || 
We show in Fig. [I] a plot of the P of Eq. || versus (3AE 
for different values of T. P becomes independent of V for 
sufficiently large magnitudes of (3AE. However, it has 
a strong T dependence around (3AE — 0, with a value 
at (3AE = of 1.0 when T = (Metropolis) and 0.0 
as r — > oo. Also, since P is a monotonically decreasing 
function of T for every value of f3AE, average equilibrated 
acceptance rates decrease monotonically with T. 

Two similar generalizations which interpolate between 
the Metropolis and symmetric decision rules were pro- 
posed recently in the context of Monte Carlo dynamics 
by Mariz, Nobre, and Tsallis j$2|. We note that the 
second of these generalizations (the so-called "geomet- 
ric unification" ) can also be further extended beyond the 
symmetric limit, reducing P when (3AE — below the 
0.5 symmetric value. We would expect the results which 



Auxiliary Field Quantum Monte Carlo 
Considerations 



The previously-cited results regarding long-run statis- 
tical errors |^],^8j can be extended for the P of Eq. || to 
any T, and would suggest that the Metropolis decision 
rule (r = 0) is usually the most efficient in that con- 
text. However, as mentioned in the Introduction, there 
is the additional consideration for auxiliary field quan- 
tum Monte Carlo that accepting a move is much 
more computationally expensive than rejecting one, with 
move acceptance dominating the computation. Then, for 
greatest efficiency, one wishes to minimize not the warm- 
up and equilibrated autocorrelation times t w and r e but 
rather the quantities A w t w and A e T e , where A w is an av- 
erage move acceptance rate during equilibration and A e 
is the average acceptance rate after equilibration. Even 
assuming that T = (Metropolis) had the smallest r's, 
it is possible that an increase in t w or r e could be over- 
compensated by a decrease in A w or A e , so that the op- 
timal r would assume some nonzero value as opposed to 
T = 0. 

Specifically, we will use the Blankenbecler-Scalapino- 
Sugar (BSS) quantum Monte Carlo algorithm to 
explore the issue of optimal T with the two-dimensional 
Hubbard model as a "test Hamiltonian" . The Hub- 
bard model interactions are decoupled using the discrete 
Ising Hubbard-Stratonovich transformation introduced 
by Hirsch ||. The Monte Carlo sampling is then per- 
formed on the Ising variables, with effective Boltzmann 
weights which depend on determinants of dense matrices 
resulting from the integrating out of the fermion degrees 
of freedom. We note that the usual Monte Carlo proce- 
dure (and the one which we will follow) is to systemat- 
ically move through the "Ising lattice" when proposing 
moves, to which the Pcskun analysis docs not rigorously 
apply |2q ]. Also, the "step size" for this type of simula- 
tion is fixed, as the only possible move is a spin flip where 
the Ising spin changes sign. Hence, we do not consider 
in this paper effects of varying step size. 

A potential additional factor is the cost required for nu- 
merically stabilizing the BSS algorithm p, ^0| , ^3[ , which 
for a given statistical error is proportional to t w and r e 
rather than A w t w and A e r e . At relatively high and inter- 
mediate temperatures this cost is negligible, but it may 
constitute a significant fraction of the computational time 
at very low temperatures. In that case, one wishes rather 
to minimize quantities of the general form t w (A w +k) and 
T e (A e + k), where the A's and k may be of comparable 
size. In the simulations which we describe, stabilization 
required a small fraction of the computational time. 



III. RESULTS 



we will discuss for the P of Eq. E 
to the generalizations of Ref. |32[ 



to apply qualitatively 
as well. 
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A. Procedure 

Our goal was to explore the behavior in auxiliary field 
quantum Monte Carlo simulations of the optimal T in 
the Monte Carlo decision rule of Eq. ||. This optimal T 
was defined as that which led to the greatest simulation 
efficiency: i.e., lowest statistical error per computational 
cost. As noted previously, we estimated this efficiency 
by the quantities A w t w during equilibration and A e r e 



J 




after equilibration, since for large systems and typical 
temperatures the leading contribution to the computa- 
tional cost scales linearly with these quantities. Again, 
A w and A w are average acceptance rates during and after 
equilibration, respectively, and t w is the equilibration or 
"warm-up" time and T e the autocorrelation time. 

Specifically, we simulated the two-dimensional Hub- 
bard model |54j , given by the Hamiltonian 



+ U ^2(n itj ,i - -)(n iii4 - -) - tx^2,ni,j, s . (7) 



Here c- , s creates an electron of spin s at site of a 

two-dimensional lattice, Jii,j, s = c\ j s Ci ; j tS is the electron 
occupation number for spin s at site and the hop- 

ping is between nearest neighbor sites. We chose hopping 
parameter t = f.O, inverse temperature (3t — t/(fcgT) = 
8, Coulomb repulsions U/t = 4 and U/t = 8, and chemi- 
cal potential ji/t = —1.2 (giving an electron density per 
site of approximately 0.76, with /i = 0.0 corresponding 
to half filling). The imaginary time step was At = 0.125, 
so that each Monte Carlo sweep contained 64 "imaginary 
time" slices. All simulations were performed on a 6 x 6 
lattice with periodic boundary conditions. 

All of the A w and t w and most of the A e and r e data 
were obtained on the massively parallel Intel Paragon at 
the Sandia National Lab MPCRL, and most of the runs 
were on the full 1824 processors. Some of the A e and 
r e data were also obtained using the CRAY T90 at the 
SDSC. 

To obtain A w and t w data, each utilized node of the 
Paragon was given a different random initial Ising config- 
uration. The acceptance rates were then collected from 
each node after each sweep through the entire (space)- 
(imaginary time) Ising lattice. The averaging of the data 
from all the nodes increased the signal-to-noise ratio to 
the extent that fits could then be made to the decay of 
the acceptance rate from its initial higher value Aq, when 
moves are more likely due to the high probability of being 
in a low- weight state, to the lower equilibrated value A e . 
A w was calculated as the average of the initial acceptance 
rate Aq and the equilibrated rate A e . Several observables 
O (n = n j + rij , a z = n j — , the staggered spin struc- 
ture functions S zz (tt,tt) and S xx (ir,w), and kinetic and 
total energies f|,|35|]) were monitored similarly to the ac- 
ceptance rate, and it was found that they all equilibrated 
at least as quickly as did the acceptance rate itself. Fits 
of the form A e + (Aq — A e ) exp(—r / t w ) were made to 
the acceptance rate data, defining the t w 's. Typically, 
only data up to twenty sweeps was used in this fit, as 
the signal became lost in the noise for a larger number of 
sweeps. 



A somewhat analogous procedure was followed in com- 
puting the r e 's. First, for the observables O listed above, 
autocorrelation functions C(t) were calculated from the 
equilibrated data using the formula 

_ (OjOj+r) 

where r is the lag time. The autocorrelations were then 
fit to the form exp(~ t/t c ), which we here took to de- 
fine the equilibrated r e autocorrelation times . Better 
statistics could be obtained for the T e 's due to the exis- 
tence of more equilibrated than equilibrating (warm-up) 
data. 

Of those considered, it was found that the observable 
with the longest calculated autocorrelation time, r e , was 
n-f— ni- Hence, we will focus primarily on the r e behavior 
of n T -nj.. 



B. Warm-Up Results 



In Fig. |2| we show the equilibration of the acceptance 
rate during warm-up for r = (Metropolis decision) and 
U/t = 4 using data averaged from 1824 processors. We 
also show an exponential decay fit to the relaxation of the 
acceptance rate from its initial to its equilibrated value. 

Such fits for various values of T were used to obtain 
the A w t w results for U/t = 4 of Fig. ||, where again 
A w = 1/2 (Aq + A e ) is an average warm-up acceptance 
rate. (The error bars for U/t — 8 were too large to 
observe statistically significant differences.) Note the 
drop shown in Fig. || in A w t w when going from T = 
(Metropolis decision rule) to T = 1 ("symmetric" rule), 
suggesting increased efficiency. 
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FIG. 2. The acceptance rate for r = and U/t = 4, shown 
starting from random configurations at step (one sweep con- 
tains 64 imaginary time steps) until it equilibrates. For clar- 
ity, error bars are shown only once every 64 points, but all 
error bars are of the same general magnitude. The points are 
averaged data from 1824 processors. The dashed line shows 
a fit of the function 0.675 + 0.123 exp(-r/l. 42) 
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FIG. 3. A w r w versus T for U/t = 4. 



C. Equilibrated Results 

In Fig. | and Fig. |, we show data for the autocorrela- 
tion C (t) of Eq. H with operator O = rtf — . The data 
were obtained by taking samples of n-\ — at every time 
slice for a large number of Monte Carlo sweeps through 
the lattice. Also given are Monte Carlo autocorrelation 
times, r e , obtained by an exponential fit, as well as the 



corresponding error bars. The r e 's and errors were cal- 
culated with the use of a Fourier transform method Q . 
As might be expected from the lower acceptance rates for 
larger F, the r e autocorrelation times increase for larger 

r. 

U/t=4 pt=8 6-by-6 |x=-1 .2 
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FIG. 4. Autocorrelation of n\ — ni for various values of V 
for U/t = 4. These curves were fit to the function exp(— r/r e ) 
to obtain r e . 
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FIG. 5. Autocorrelation of n-f — for various values of V 
for U/t = 8. These curves were fit to the function exp(— r/r e ) 
to obtain r e . 

As shown in Fig ^, we observed a "ringing" in the au- 
tocorrelation for the staggered spin structure function in 
the z-direction, S zz (ir,ir). No such ringing was observed 
for any other observable studied, including the staggered 
spin structure function in the x-direction, S xx (tt,tt). (It 
is known that a BSS simulation with discrete Hubbard- 
Stratonovich decoupling can lead to different variances 
for quantities with the same averages which are measured 
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in different directions p5|.) 



IV. DISCUSSION 



U/t=4 (3t=8 6-by-6 n=-1.2 



1.0 



r=o x e =o.3±o.i 

T=2 x e =1.4±0.3 

r=8 x =4.1 ± 0.4 



We have explored possible improvements in the effi- 
ciency of determinantal fermion Monte Carlo simulations 
using a generalization of the standard Metropolis decision 
rule, 
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FIG. 6. Autocorrelation of S zz (n,n) for various values of 
r for U/t — 4. Taking just the peaks, these curves can also 
be fit to exp(— r/r e ) to obtain a r e for this observable. 
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FIG. 7. A e r e versus F for J7/t = 4 and U/t = 8. 

In Fig. [?], we show A e r e for both U/t — 4 (top) and 
U/t = 8 (bottom) . There were no statistically significant 
changes observed in A e r e for U/t = 8 in the studied range 
< r < 5. However, for U/t = 4, A e r e was observed to 
increase with increasing T, and a monotonic increase is 
quite consistent with the data. Of the < T < 8 values 
studied for U/t = 4, V — gave the (statistically signifi- 
cant) lowest value of A e r e , suggesting that the Metropolis 
decision was most efficient in this case. 
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(JAE 



l+re 



r+e- 



AE > 
— , AE < . 



(9) 



This generalization interpolates between the Metropo- 
lis |I| and so-called "symmetric" |^,^6j decision rules 
but also has additional flexibility: the average accep- 
tance rate can be smoothly reduced from the (maximal) 
Metropolis value (when T = 0) through the symmet- 
ric value (r = 1) down to an arbitrarily small value. 
Specifically, we performed auxiliary field quantum Monte 
Carlo simulations [p|-|lo|; however, our approach and re- 
sults may have application to electronic structure fermion 
Monte Carlo methods as well |11-16|. 

One might in general expect that a lower acceptance 
rate would be less computationally efficient, since the 
sampling would then proceed more slowly through the 
configuration space. However, determinantal quantum 
Monte Carlo simulations have the feature that an ac- 
cepted move is generally much more costly than a re- 
jected one and that the calculations associated with 
move acceptance will typically dominate the computa- 
tion. Then, instead of trying to minimize autocorrela- 
tion times r e for greater efficiency, one rather wishes to 
minimize A e r e , where A e is the average equilibrated ac- 
ceptance rate. A similar argument holds for A w t w as op- 
posed to t w , where t w is the equilibration or "warm-up" 
time and A w is the averaged acceptance rate during the 
equilibration process. Since equilibration poses a poten- 
tial (serial) bottleneck in parallelization, we considered 
the behavior both of A e r e and of A w t w . In particular, 
we explored whether generalizations of the Metropolis 
decision could reduce A e r e or A w r w . 

For a test system, we utilized the two-dimensional 
Hubbard model at approximately 3/4 filling with U/t = 4 
(weak coupling) and U/t = 8 (moderate to strong cou- 
pling). For computing an autocorrelation time r e , we 
used the observable n-\ — nj_, which had the longest such 
calculated autocorrelation time of any of the observables 
considered. For calculating the "warm-up" time t w , we 
monitored how the acceptance rate equilibrated from ran- 
dom initial configurations. When U/t — 4, both r e and 
A e r e increased with T, suggesting that T = (Metropolis 
decision) was the most efficient. The r e 's grew monotoni- 
cally with T when U/t = 8, but there was no statistically 
significant variation observed in A e r e for the T's tested. 
The error bars in A w t w when U/t = 8 were too large for 
meaningful comparisons; however, when U/t = 4, A w t w 
was approximately halved in going from T — (Metropo- 
lis) to r = 1 (symmetric), indicating that the symmetric 
decision was more efficient during equilibration. 
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A byproduct of this work was calculations of autocorre- 
lation times r e and "warm-up" times t w for some sample 
auxiliary field (two-dimensional Hubbard model) quan- 
tum Monte Carlo runs, which calculations to the best of 
our knowledge have not previously appeared. For both 
Metropolis and symmetric decision rules, the r's were 
typically at most a few sweeps through the Ising (space)- 
(imaginary-time) auxiliary field lattice. These values are 
shorter than might have been expected, and suggest that 
the often several thousand warm-up sweeps convention- 
ally done in vector simulations are very adequate. Such 
r values may also provide useful "ball park" estimates 
when planning parallel runs. 

It was clear from our simulations that which of the 
Monte Carlo decisions is most efficient in determinantal 
fcrmion Monte Carlo can depend on specifics of the model 
and parameters. However, the standard Metropolis deci- 
sion was optimal or near optimal once equilibration had 
been reached in the two sample cases studied, suggest- 
ing that it is efficient after warm-up. The symmetric 
decision, however, was more efficient during warm-up. 
Particularly if the (serial) warm-up process becomes a 
bottleneck in parallel simulations, this suggests that the 
use of symmetric or other non-Metropolis decision rules 
could lead to greater computational speed. 
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